traits
above- and belowground
acquisitive vs. conservative resource-use strategies
collaborative vs. individual resource-use strategies
fine root proportion
thickness of roots
Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.
Author: Caryn D. Iwanaga
Updated: 02/06/2025
library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)
cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
mutate(
nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))
# Define a reusable color scale
# custom_colors <- scale_fill_manual(
# values = c(
# "control" = rgb(195, 197, 193, maxColorValue = 250),
# "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
# "+compost" = rgb(167, 156, 109, maxColorValue = 225)
# )
# )
custom_colors <- scale_fill_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
# Set a global theme
theme_set(theme_bw())
site.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
mutate(
grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
) %>% # separate by block (grazing intensities)
group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
matrix0 <- spread(
site.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix <- matrix0[, c(5:79)]
# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)
site.block.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
summarise(
# pct_cover,
abundance = mean(pct_cover)
# sd = sd(pct_cover)
)
`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.
site.block.sd.df <-
cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
summarise(
pct_cover,
abundance = mean(pct_cover),
sd = sd(pct_cover)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.
matrix0.block <- spread(
site.block.df,
code4,
abundance,
fill = 0)
# make rownames ----
rownames(matrix0.block) <- paste(matrix0.block$nut_trt, matrix0.block$ppt_trt, matrix0.block$yr, matrix0.block$block, sep = "_")
Warning: Setting row names on a tibble is deprecated.
# delete columns ----
matrix.block <- matrix0.block[, -c(0:4)]
# relativize data with Bray-Curtis dissimilarity
matrix.block.bray <- decostand(matrix.block, "total")
rownames(matrix.block.bray) <- rownames(matrix0.block)
Rules of thumb:
< 0.05 is excellent
0.05-0.1 is good
0.1-0.2 is fair
0.2-0.3 is cause for concern…
> 0.3 is poor, random
Stress: fair
0.1966553
Grouping by year and grazing history had most significant clustering shown.
nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)
Run 0 stress 0.2020736
Run 1 stress 0.2037794
Run 2 stress 0.1973259
... New best solution
... Procrustes: rmse 0.07504543 max resid 0.2323388
Run 3 stress 0.1970656
... New best solution
... Procrustes: rmse 0.02167554 max resid 0.1184388
Run 4 stress 0.2112285
Run 5 stress 0.2016301
Run 6 stress 0.2034469
Run 7 stress 0.2100941
Run 8 stress 0.2020249
Run 9 stress 0.1990567
Run 10 stress 0.205915
Run 11 stress 0.2072425
Run 12 stress 0.2077976
Run 13 stress 0.196667
... New best solution
... Procrustes: rmse 0.01874332 max resid 0.1106696
Run 14 stress 0.2021165
Run 15 stress 0.2061482
Run 16 stress 0.1966667
... New best solution
... Procrustes: rmse 0.0005495779 max resid 0.002877469
... Similar to previous best
Run 17 stress 0.2158361
Run 18 stress 0.21463
Run 19 stress 0.2120927
Run 20 stress 0.2037336
*** Best solution repeated 1 times
nmds.comp # pretty not great stress
Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray
Distance: bray
Dimensions: 2
Stress: 0.1966667
Stress type 1, weak ties
Best solution was repeated 1 time in 20 tries
The best solution was from try 16 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray’
stressplot(nmds.comp)
plot(nmds.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))
nmds_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: cause for concern
0.2329306
nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)
Run 0 stress 0.2329306
Run 1 stress 0.2385405
Run 2 stress 0.2608291
Run 3 stress 0.2618978
Run 4 stress 0.2625943
Run 5 stress 0.2369671
Run 6 stress 0.2329931
... Procrustes: rmse 0.002596777 max resid 0.01287822
Run 7 stress 0.2447034
Run 8 stress 0.2406908
Run 9 stress 0.2476092
Run 10 stress 0.2412944
Run 11 stress 0.232978
... Procrustes: rmse 0.002270404 max resid 0.012997
Run 12 stress 0.23922
Run 13 stress 0.2497086
Run 14 stress 0.2478848
Run 15 stress 0.2365873
Run 16 stress 0.2329684
... Procrustes: rmse 0.001740474 max resid 0.01298233
Run 17 stress 0.2520339
Run 18 stress 0.248187
Run 19 stress 0.2380139
Run 20 stress 0.2333533
... Procrustes: rmse 0.0115742 max resid 0.112582
Run 21 stress 0.2489426
Run 22 stress 0.2400835
Run 23 stress 0.2421075
Run 24 stress 0.2510936
Run 25 stress 0.232944
... Procrustes: rmse 0.001509127 max resid 0.01190131
*** Best solution was not repeated -- monoMDS stopping criteria:
24: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
nmds.block.comp
Call:
metaMDS(comm = matrix.block.bray, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.block.bray
Distance: bray
Dimensions: 2
Stress: 0.2329306
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.block.bray’
stressplot(nmds.block.comp)
plot(nmds.block.comp, type="t")
nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))
nmds_dat.block <- cbind(nmds1, nmds2) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year - Blocks")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment - Blocks") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(block), color = as.factor(block))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Block")
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(block), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(block)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Block & Year - Blocks") +
scale_color_manual(
values = c(
# "low" = "cyan",
# "high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
Stress: fair
0.1869068
# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
filter(grazing_hist == "low")
row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")
matrix.low <- matrix.low[, c(5:79)]
rownames(matrix.low) <- row.low
Warning: Setting row names on a tibble is deprecated.
matrix.bray.low <- decostand(matrix.low, "total")
nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)
Run 0 stress 0.2000574
Run 1 stress 0.1960027
... New best solution
... Procrustes: rmse 0.03759248 max resid 0.136796
Run 2 stress 0.2299212
Run 3 stress 0.1908791
... New best solution
... Procrustes: rmse 0.1110161 max resid 0.2532458
Run 4 stress 0.1869068
... New best solution
... Procrustes: rmse 0.1544638 max resid 0.4035849
Run 5 stress 0.2066946
Run 6 stress 0.1878648
Run 7 stress 0.2004585
Run 8 stress 0.1876855
Run 9 stress 0.2113763
Run 10 stress 0.2243075
Run 11 stress 0.1964757
Run 12 stress 0.1893384
Run 13 stress 0.2100567
Run 14 stress 0.2455646
Run 15 stress 0.1930573
Run 16 stress 0.2043372
Run 17 stress 0.245567
Run 18 stress 0.2240266
Run 19 stress 0.1942011
Run 20 stress 0.1893384
Run 21 stress 0.2104986
Run 22 stress 0.1960027
Run 23 stress 0.22813
Run 24 stress 0.2593359
Run 25 stress 0.1908791
*** Best solution was not repeated -- monoMDS stopping criteria:
23: stress ratio > sratmax
2: scale factor of the gradient < sfgrmin
nmds.comp.low # pretty not great stress
Call:
metaMDS(comm = matrix.bray.low, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.low
Distance: bray
Dimensions: 2
Stress: 0.1869068
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 4 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.low’
stressplot(nmds.comp.low)
plot(nmds.comp.low, type="t")
nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
)
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))
nmds_low_dat <- cbind(nmds1, nmds2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, color = as.factor(yr), shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# add year ellipses ----
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_low_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
custom_colors +
labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
filter(yr == 2019)
row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")
matrix.2019 <- matrix.2019[, c(5:79)]
rownames(matrix.2019) <- row.2019
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2019 <- decostand(matrix.2019, "total")
Stress: Excellent
0.0560178
nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)
Run 0 stress 0.0560178
Run 1 stress 0.06223133
Run 2 stress 0.05984339
Run 3 stress 0.05905532
Run 4 stress 0.06384512
Run 5 stress 0.05851802
Run 6 stress 0.05852022
Run 7 stress 0.06942148
Run 8 stress 0.05654879
Run 9 stress 0.06942148
Run 10 stress 0.05654879
Run 11 stress 0.05877078
Run 12 stress 0.0662005
Run 13 stress 0.05982901
Run 14 stress 0.05957995
Run 15 stress 0.07018407
Run 16 stress 0.05652852
Run 17 stress 0.06074817
Run 18 stress 0.05652852
Run 19 stress 0.05654879
Run 20 stress 0.07335647
Run 21 stress 0.05851807
Run 22 stress 0.05652852
Run 23 stress 0.0735608
Run 24 stress 0.0761203
Run 25 stress 0.06228134
*** Best solution was not repeated -- monoMDS stopping criteria:
13: stress ratio > sratmax
12: scale factor of the gradient < sfgrmin
nmds.2019 # pretty not great stress
Call:
metaMDS(comm = matrix.bray.2019, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2019
Distance: bray
Dimensions: 2
Stress: 0.0560178
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2019’
stressplot(nmds.2019)
plot(nmds.2019, type="t")
nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
)
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))
nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 3, stroke = 1) +
labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
labs(title = "2019 - Grouping by Nutrient Treatment")+
custom_colors
ggplot(nmds_2019_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2019 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
filter(yr == 2020)
row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")
matrix.2020 <- matrix.2020[, c(5:79)]
rownames(matrix.2020) <- row.2020
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2020 <- decostand(matrix.2020, "total")
Stress: fair
0.1522495
nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)
Run 0 stress 0.1524241
Run 1 stress 0.1584809
Run 2 stress 0.1584809
Run 3 stress 0.1688381
Run 4 stress 0.1625903
Run 5 stress 0.153917
Run 6 stress 0.3515938
Run 7 stress 0.1692911
Run 8 stress 0.1692911
Run 9 stress 0.1584809
Run 10 stress 0.1522495
... New best solution
... Procrustes: rmse 0.02107749 max resid 0.06346372
Run 11 stress 0.1584809
Run 12 stress 0.1522045
... New best solution
... Procrustes: rmse 0.01183105 max resid 0.0373035
Run 13 stress 0.1744238
Run 14 stress 0.1524241
... Procrustes: rmse 0.02481397 max resid 0.06599247
Run 15 stress 0.1742461
Run 16 stress 0.1595062
Run 17 stress 0.1524241
... Procrustes: rmse 0.02481007 max resid 0.06599555
Run 18 stress 0.1625903
Run 19 stress 0.1522495
... Procrustes: rmse 0.01183378 max resid 0.03737577
Run 20 stress 0.1522495
... Procrustes: rmse 0.0118299 max resid 0.03737025
Run 21 stress 0.1747777
Run 22 stress 0.184858
Run 23 stress 0.1524241
... Procrustes: rmse 0.02481216 max resid 0.06599381
Run 24 stress 0.1866314
Run 25 stress 0.2912948
*** Best solution was not repeated -- monoMDS stopping criteria:
20: stress ratio > sratmax
5: scale factor of the gradient < sfgrmin
nmds.2020
Call:
metaMDS(comm = matrix.bray.2020, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2020
Distance: bray
Dimensions: 2
Stress: 0.1522045
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 12 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2020’
stressplot(nmds.2020)
plot(nmds.2020, type="t")
nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
)
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))
nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 1) +
labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Nutrient Treatment")
ggplot(nmds_2020_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2020 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
filter(yr == 2021)
row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")
matrix.2021 <- matrix.2021[, c(5:79)]
rownames(matrix.2021) <- row.2021
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021 <- decostand(matrix.2021, "total")
Stress
0.08067789
# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "low")
row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")
matrix.2021.1 <- matrix.2021.1[, c(5:79)]
rownames(matrix.2021.1) <- row.2021.1
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")
nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)
Run 0 stress 0.08067789
Run 1 stress 0.1407365
Run 2 stress 0.131677
Run 3 stress 0.08067789
... New best solution
... Procrustes: rmse 4.448434e-06 max resid 7.224535e-06
... Similar to previous best
Run 4 stress 0.1393679
Run 5 stress 0.08067789
... New best solution
... Procrustes: rmse 1.053765e-06 max resid 1.585547e-06
... Similar to previous best
Run 6 stress 0.1393365
Run 7 stress 0.1393365
Run 8 stress 0.08067789
... Procrustes: rmse 5.054553e-06 max resid 8.236855e-06
... Similar to previous best
Run 9 stress 0.08067789
... Procrustes: rmse 5.780641e-06 max resid 8.023545e-06
... Similar to previous best
Run 10 stress 0.08067789
... Procrustes: rmse 4.42502e-06 max resid 6.754172e-06
... Similar to previous best
Run 11 stress 0.2039408
Run 12 stress 0.08067789
... Procrustes: rmse 1.729526e-06 max resid 2.63632e-06
... Similar to previous best
Run 13 stress 0.08067789
... New best solution
... Procrustes: rmse 1.098472e-06 max resid 1.772254e-06
... Similar to previous best
Run 14 stress 0.1316771
Run 15 stress 0.2639821
Run 16 stress 0.203941
Run 17 stress 0.08067789
... Procrustes: rmse 3.956956e-06 max resid 5.336733e-06
... Similar to previous best
Run 18 stress 0.08067789
... Procrustes: rmse 1.630904e-06 max resid 3.249717e-06
... Similar to previous best
Run 19 stress 0.1203381
Run 20 stress 0.2039414
*** Best solution repeated 3 times
nmds.2021.1
Call:
metaMDS(comm = matrix.bray.2021.1, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021.1
Distance: bray
Dimensions: 2
Stress: 0.08067789
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 13 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021.1’
stressplot(nmds.2021.1)
plot(nmds.2021.1, type="t")
dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))
nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))
nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
# geom_point(size = 4, stroke = 0.5) +
# labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")
Stress
0.08065034
# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
filter(yr == 2021,
grazing_hist == "high")
row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")
matrix.2021.2 <- matrix.2021.2[, c(5:79)]
rownames(matrix.2021.2) <- row.2021.2
Warning: Setting row names on a tibble is deprecated.
matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")
nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)
Run 0 stress 0.08065034
Run 1 stress 0.1378308
Run 2 stress 0.08065034
... Procrustes: rmse 2.805412e-05 max resid 5.594612e-05
... Similar to previous best
Run 3 stress 0.1378308
Run 4 stress 0.08065035
... Procrustes: rmse 7.802029e-05 max resid 0.0001548812
... Similar to previous best
Run 5 stress 0.08065034
... Procrustes: rmse 3.731121e-05 max resid 6.566436e-05
... Similar to previous best
Run 6 stress 0.1378308
Run 7 stress 0.08554815
Run 8 stress 0.08065037
... Procrustes: rmse 0.0001137107 max resid 0.0002575052
... Similar to previous best
Run 9 stress 0.08554815
Run 10 stress 0.1378309
Run 11 stress 0.08065034
... New best solution
... Procrustes: rmse 1.388738e-05 max resid 1.990997e-05
... Similar to previous best
Run 12 stress 0.1378309
Run 13 stress 0.1905268
Run 14 stress 0.1743257
Run 15 stress 0.08065036
... Procrustes: rmse 7.179091e-05 max resid 0.0001453003
... Similar to previous best
Run 16 stress 0.08065034
... Procrustes: rmse 4.443671e-05 max resid 7.380012e-05
... Similar to previous best
Run 17 stress 0.08065034
... Procrustes: rmse 9.112915e-06 max resid 1.854915e-05
... Similar to previous best
Run 18 stress 0.08065035
... Procrustes: rmse 4.869656e-05 max resid 9.470507e-05
... Similar to previous best
Run 19 stress 0.1460913
Run 20 stress 0.1378308
*** Best solution repeated 5 times
nmds.2021.2
Call:
metaMDS(comm = matrix.bray.2021.2, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021.2
Distance: bray
Dimensions: 2
Stress: 0.08065034
Stress type 1, weak ties
Best solution was repeated 5 times in 20 tries
The best solution was from try 11 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021.2’
stressplot(nmds.2021.2)
plot(nmds.2021.2, type="t")
dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))
nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
mutate(trt = rownames(dummy)) %>%
separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")
nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))
nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")
Stress: fair
0.1230831
nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)
Run 0 stress 0.1232098
Run 1 stress 0.1339879
Run 2 stress 0.1230831
... New best solution
... Procrustes: rmse 0.04463789 max resid 0.1321469
Run 3 stress 0.1230831
... Procrustes: rmse 1.350077e-05 max resid 4.219102e-05
... Similar to previous best
Run 4 stress 0.133988
Run 5 stress 0.1507047
Run 6 stress 0.1573377
Run 7 stress 0.1230831
... Procrustes: rmse 2.810056e-06 max resid 8.291934e-06
... Similar to previous best
Run 8 stress 0.1339879
Run 9 stress 0.1230831
... Procrustes: rmse 8.646086e-06 max resid 2.745517e-05
... Similar to previous best
Run 10 stress 0.132972
Run 11 stress 0.1232098
... Procrustes: rmse 0.04463824 max resid 0.1315186
Run 12 stress 0.1230831
... New best solution
... Procrustes: rmse 2.871318e-06 max resid 7.467703e-06
... Similar to previous best
Run 13 stress 0.1421227
Run 14 stress 0.1232098
... Procrustes: rmse 0.04463976 max resid 0.1315202
Run 15 stress 0.1339879
Run 16 stress 0.1515614
Run 17 stress 0.1230832
... Procrustes: rmse 2.468942e-05 max resid 7.756011e-05
... Similar to previous best
Run 18 stress 0.1573377
Run 19 stress 0.1230831
... Procrustes: rmse 7.42768e-06 max resid 2.388952e-05
... Similar to previous best
Run 20 stress 0.132972
*** Best solution repeated 3 times
nmds.2021
Call:
metaMDS(comm = matrix.bray.2021, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.2021
Distance: bray
Dimensions: 2
Stress: 0.1230831
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 12 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.2021’
stressplot(nmds.2021)
plot(nmds.2021, type="t")
nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
)
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))
nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
select(nut_trt:grazing_hist, NMDS1, NMDS2)
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, color = grazing_hist, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block <- matrix0.block %>%
filter(yr == 2021)
row.2021.block <- paste(matrix.2021.block$nut_trt, matrix.2021.block$ppt_trt, matrix.2021.block$yr, matrix.2021.block$block, sep = "_")
matrix.2021.block <- matrix.2021.block[, c(5:79)]
rownames(matrix.2021.block) <- row.2021.block
Warning: Setting row names on a tibble is deprecated.
matrix.bray.block.2021 <- decostand(matrix.2021.block, "total")
Stress: fair
0.1702944
nmds.2021.block <- metaMDS(matrix.bray.block.2021, k=2, trymax = 25)
Run 0 stress 0.1702944
Run 1 stress 0.2323726
Run 2 stress 0.2333131
Run 3 stress 0.1702944
... New best solution
... Procrustes: rmse 8.12052e-06 max resid 2.925004e-05
... Similar to previous best
Run 4 stress 0.1704654
... Procrustes: rmse 0.006480506 max resid 0.02760662
Run 5 stress 0.2085283
Run 6 stress 0.2044772
Run 7 stress 0.2146073
Run 8 stress 0.1942847
Run 9 stress 0.1873962
Run 10 stress 0.2098049
Run 11 stress 0.1917675
Run 12 stress 0.1833145
Run 13 stress 0.2044772
Run 14 stress 0.226506
Run 15 stress 0.2142311
Run 16 stress 0.1702944
... Procrustes: rmse 5.781581e-06 max resid 1.825294e-05
... Similar to previous best
Run 17 stress 0.1798315
Run 18 stress 0.1704654
... Procrustes: rmse 0.006481124 max resid 0.02760841
Run 19 stress 0.2355275
Run 20 stress 0.2089454
*** Best solution repeated 2 times
nmds.2021.block
Call:
metaMDS(comm = matrix.bray.block.2021, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021
Distance: bray
Dimensions: 2
Stress: 0.1702944
Stress type 1, weak ties
Best solution was repeated 2 times in 20 tries
The best solution was from try 3 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021’
stressplot(nmds.2021.block)
plot(nmds.2021.block, type="t")
nmds1.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 4)
)
nmds2.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(2), display=c("sites")))
nmds_2021_dat.block <- cbind(nmds1.2021.block, nmds2.2021.block) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
filter(yr == 2021,
block == 3 | block == 4)
row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")
matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]
rownames(matrix.2021.block.low) <- row.2021.block.low
Warning: Setting row names on a tibble is deprecated.
matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")
Stress: fair
0.1282296
nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)
Run 0 stress 0.1282296
Run 1 stress 0.1282296
... New best solution
... Procrustes: rmse 6.900845e-05 max resid 0.0001846812
... Similar to previous best
Run 2 stress 0.1282296
... New best solution
... Procrustes: rmse 6.471065e-05 max resid 0.0002196127
... Similar to previous best
Run 3 stress 0.1282296
... Procrustes: rmse 8.428118e-05 max resid 0.0003017723
... Similar to previous best
Run 4 stress 0.1985506
Run 5 stress 0.1282296
... Procrustes: rmse 1.033532e-05 max resid 2.797048e-05
... Similar to previous best
Run 6 stress 0.1785638
Run 7 stress 0.1282296
... Procrustes: rmse 8.012815e-05 max resid 0.0002839093
... Similar to previous best
Run 8 stress 0.1282296
... Procrustes: rmse 5.434575e-05 max resid 0.0001894526
... Similar to previous best
Run 9 stress 0.1282296
... New best solution
... Procrustes: rmse 1.600545e-05 max resid 5.619322e-05
... Similar to previous best
Run 10 stress 0.1282296
... Procrustes: rmse 1.523127e-05 max resid 5.391544e-05
... Similar to previous best
Run 11 stress 0.1428149
Run 12 stress 0.1627541
Run 13 stress 0.2101121
Run 14 stress 0.1282296
... Procrustes: rmse 2.043501e-05 max resid 6.96312e-05
... Similar to previous best
Run 15 stress 0.2053314
Run 16 stress 0.2106634
Run 17 stress 0.1282296
... Procrustes: rmse 8.924757e-05 max resid 0.0003091915
... Similar to previous best
Run 18 stress 0.1428149
Run 19 stress 0.1428149
Run 20 stress 0.1428149
*** Best solution repeated 4 times
nmds.2021.block.low
Call:
metaMDS(comm = matrix.bray.block.2021.low, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.low
Distance: bray
Dimensions: 2
Stress: 0.1282296
Stress type 1, weak ties
Best solution was repeated 4 times in 20 tries
The best solution was from try 9 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.low’
stressplot(nmds.2021.block.low)
plot(nmds.2021.block.low, type="t")
nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
)
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))
nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")
# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.high <- matrix0.block %>%
filter(yr == 2021,
block == 1 | block == 2)
row.2021.block.high <- paste(matrix.2021.block.high$nut_trt, matrix.2021.block.high$ppt_trt, matrix.2021.block.high$yr, matrix.2021.block.high$block, sep = "_")
matrix.2021.block.high <- matrix.2021.block.high[, -c(0:4)]
rownames(matrix.2021.block.high) <- row.2021.block.high
Warning: Setting row names on a tibble is deprecated.
matrix.bray.block.2021.high <- decostand(matrix.2021.block.high, "total")
Stress: fair
0.1412997
nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)
Run 0 stress 0.1412997
Run 1 stress 0.1430823
Run 2 stress 0.1430824
Run 3 stress 0.1856239
Run 4 stress 0.1581575
Run 5 stress 0.1412999
... Procrustes: rmse 0.0003384719 max resid 0.001097121
... Similar to previous best
Run 6 stress 0.1430825
Run 7 stress 0.1581574
Run 8 stress 0.1856239
Run 9 stress 0.1412997
... Procrustes: rmse 3.643604e-05 max resid 0.0001217465
... Similar to previous best
Run 10 stress 0.1581574
Run 11 stress 0.1581574
Run 12 stress 0.1430824
Run 13 stress 0.1430824
Run 14 stress 0.1581574
Run 15 stress 0.2412356
Run 16 stress 0.1430827
Run 17 stress 0.1430824
Run 18 stress 0.1412998
... Procrustes: rmse 0.0002267585 max resid 0.0007337562
... Similar to previous best
Run 19 stress 0.1581575
Run 20 stress 0.1581574
*** Best solution repeated 3 times
nmds.2021.block.high
Call:
metaMDS(comm = matrix.bray.block.2021.high, k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray.block.2021.high
Distance: bray
Dimensions: 2
Stress: 0.1412997
Stress type 1, weak ties
Best solution was repeated 3 times in 20 tries
The best solution was from try 0 (metric scaling or null solution)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray.block.2021.high’
stressplot(nmds.2021.block.high)
plot(nmds.2021.block.high, type="t")
nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
mutate(
nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
)
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))
nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
select(nut_trt:block, NMDS1, NMDS2)
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, color = block, shape = nut_trt)) +
geom_point(size = 4, stroke = 0.5) +
labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")
ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2, fill = ppt_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
custom_colors +
labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")
Greenhouse traits
Field traits
“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”
“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”
Workflow plan:
overlay species trait data as vectors on the NMDS
only greenhouse traits available (only sla data from field)
use PC scores as composite trait axes OR
average the traits for each species
# # read in greenhouse trait data
# traits.gh <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=ya1ck6p7&dl=1")
#
# # read in field trait data
# # download.file("https://www.dropbox.com/scl/fi/trqrnz54dkjoa9do4tle5/leaf-area-data_NG.xlsx?rlkey=s0kic2a5yqknikc4icek8y7rn&st=dl4x634h&dl=1",destfile = "field_sla.xlsx")
# traits.sla.field <-
# read_excel("field_sla.xlsx",col_names = TRUE,col_types = "text")
#
# rda(traits.gh, scale = TRUE)
Trait data is identical. Will clean data in
# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
#
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
#
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)
traits.cleaned <- read.csv("data/GreenhouseTraits_corrected_20240212.csv") %>%
mutate(
Fresh.leaf.mass..g. = as.numeric(ifelse(
Fresh.leaf.mass..g. == "Skipped accidentally", NA, Fresh.leaf.mass..g.
)),
Dry.leaf.mass..g. = as.numeric(ifelse(
Dry.leaf.mass..g. == "No sample", NA, Dry.leaf.mass..g.
)),
Root.dry.biomass..g. = as.numeric(ifelse(
Root.dry.biomass..g. == "roots lost", NA, Root.dry.biomass..g.
)),
Root.volume..cm3. = as.numeric(ifelse(
Root.volume..cm3. == "not scanned", NA, Root.volume..cm3.
))
)
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `Root.dry.biomass..g. = as.numeric(ifelse(Root.dry.biomass..g. == "roots lost", NA,
Root.dry.biomass..g.))`.
Caused by warning:
! NAs introduced by coercion
# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
mutate(
ID = ifelse(
ID == "Agoseris", "AGSP", ifelse(
ID == "AVEBAR", "AVBA", ifelse(
ID == "BRDR", "BRDI", ID
)
)
)
)
# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
group_by(ID) %>%
summarise(
across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
)))
# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)
intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)
[1] "ASSP" "AST1" "AST2" "AST3" "CIQU" "CRTI" "CYEC" "DACA" "DIMU" "GAGR" "GAPA" "GAPH" "HYSP"
[14] "JUBU" "LYHY" "MAD1" "MAD2" "NAPU" "PAVI" "SABI1" "SABI2" "SIGA" "THGR" "TRSP" "UNBU" "UNF1"
[27] "UNF3" "UNF4" "UNF8" "UNFR" "UNGR" "XAST"
setdiff(traits.means$ID, colnames(matrix.bray))
[1] "ACMI" "ANARp" "AVFA" "BRCA" "BRHOp" "BRNIf" "BRNIp" "CLPUf" "CLPUp"
[10] "CYDA" "Crassula" "ESCA" "FEMI" "FEMY" "GITRf" "GITRp" "HYGL" "LACA"
[19] "LENIf" "LENIp" "LOPU" "LUBI" "Leontodon" "MAELf" "MAELp" "MICA" "NEMA"
[28] "PLERf" "PLERp" "RUPU" "Sanicula" "TRCI" "TRHIpi" "TRWIf" "TRWIp" "TRWIpi"
[37] "VIVA"
# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
dplyr::select(code4, species) %>%
summarise(
code4 = unique(code4),
species = unique(species)
)
Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.
unique.trait.species <- traits.cleaned %>%
dplyr::select(Taxon, ID) %>%
distinct(ID, .keep_all = TRUE)
change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)
What’s the difference between:
‘ANAR’ and ‘ANARp’ in trait data?
‘BRHO’ and ‘BRHOp’ in trait data?
‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)
‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)
‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)
‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)
‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)
‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)
‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)
‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)
‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)
‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)
‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)
change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)
change BRDR –> BRDI
Stress:
0.2011127
# subset nmds to match species names
nmds.comp.trait <- metaMDS(matrix.bray %>%
dplyr::select(intersect.names),
k=2, trymax = 25)
Run 0 stress 0.2060005
Run 1 stress 0.2007525
... New best solution
... Procrustes: rmse 0.06486829 max resid 0.1521338
Run 2 stress 0.2133228
Run 3 stress 0.2079153
Run 4 stress 0.2051687
Run 5 stress 0.2107849
Run 6 stress 0.2078775
Run 7 stress 0.2027188
Run 8 stress 0.2013372
Run 9 stress 0.2050098
Run 10 stress 0.2017137
Run 11 stress 0.2042866
Run 12 stress 0.2051682
Run 13 stress 0.2185014
Run 14 stress 0.2078775
Run 15 stress 0.2156872
Run 16 stress 0.2221242
Run 17 stress 0.2079153
Run 18 stress 0.2023083
Run 19 stress 0.2181406
Run 20 stress 0.2010909
... Procrustes: rmse 0.08100757 max resid 0.2166622
Run 21 stress 0.2086475
Run 22 stress 0.2046958
Run 23 stress 0.2061378
Run 24 stress 0.2020872
Run 25 stress 0.2080994
*** Best solution was not repeated -- monoMDS stopping criteria:
24: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
nmds.comp.trait # pretty not great stress
Call:
metaMDS(comm = matrix.bray %>% dplyr::select(intersect.names), k = 2, trymax = 25)
global Multidimensional Scaling using monoMDS
Data: matrix.bray %>% dplyr::select(intersect.names)
Distance: bray
Dimensions: 2
Stress: 0.2007525
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 1 (random start)
Scaling: centring, PC rotation, halfchange scaling
Species: expanded scores based on ‘matrix.bray %>% dplyr::select(intersect.names)’
stressplot(nmds.comp.trait)
plot(nmds.comp.trait, type="t")
nmds1 <- as.data.frame(scores(nmds.comp.trait, choices=c(1), display=c("sites")))
nmds2 <- as.data.frame(scores(nmds.comp.trait, choices=c(2), display=c("sites")))
nmds_dat_trait <- cbind(nmds1, nmds2) %>%
as.data.frame() %>% # Ensure it's a data frame
rownames_to_column(var = "site_id") %>% # Move row names into a new column
separate(site_id, into = c("nut_trt", "ppt_trt", "yr", "grazing_hist"), sep = "_")
# subset matrix and trait data to match names
traits.means.comp <- traits.means %>%
filter(ID %in% intersect.names) %>%
rename(
Height = Height..cm._mean,
FreshLeafMass = Fresh.leaf.mass..g._mean,
DryLeafMass = Dry.leaf.mass..g._mean,
LDMC = LDMC_mean, # Leaf dry matter content (LDMC, the ratio of leaf dry mass to fresh mass)
LeafArea = Leaf.Area..cm2._mean,
SLA = SLA..cm2.g._mean,
ShootDryBiomass = Shoot.dry.biomass..g._mean,
RootDryBiomass = Root.dry.biomass..g._mean,
TotalBiomass = Total.biomass..g._mean,
RMF = RMF_mean, # Root mass fraction (RMF) is a plant trait that measures the proportion of a plant's dry mass that is in its roots
RootVol = Root.volume..cm3._mean,
RootDensity = Root.density..g.cm3._mean,
CoarseRootDiameter = Coarse.root.diameter..mm._mean,
RootLength = Length..mm._mean,
FineRootLength = Fine.root.length..mm._mean,
CoarseRootLength = Coarse.root.length..mm._mean,
CoarseRootSpecLength = Coarse.root.specific.length..cm.g._mean,
FineRootSpecLength = Fine.root.specific.length..cm.g._mean,
PropFineRoots = Proportion.fine.roots_mean
)
# traits.means.comp2 <- traits.means.comp %>%
# dplyr::select(-ID)
rownames(traits.means.comp2) <- traits.means.comp$ID
Warning: Setting row names on a tibble is deprecated.
# trait vectors onto the NMDS
# envfit(nmds.comp.trait, traits.means.comp, permutations = 999)
# example from stack overflow: envfit(ord ~ var1 + var2 + var3, dune.spec, display="sp"), 'ord' is the nmds model, 'var1' is the trait1 i think?
colnames(traits.means.comp)
[1] "ID" "Height" "FreshLeafMass" "DryLeafMass"
[5] "LDMC" "LeafArea" "SLA" "ShootDryBiomass"
[9] "RootDryBiomass" "TotalBiomass" "RMF" "RootVol"
[13] "RootDensity" "CoarseRootDiameter" "RootLength" "FineRootLength"
[17] "CoarseRootLength" "CoarseRootSpecLength" "FineRootSpecLength" "PropFineRoots"
vectors.nmds.comp.trait <- envfit(nmds.comp.trait ~ Height +
LeafArea +
RMF, traits.means.comp, display = "sp") # ONLY Height, Leaf area, and RMF significant
# is it different when you only do certain traits (above- vs. below-ground)? or can you just throw everything in there and it will tell you what is significant and you can subset from there?
# Looks like different subsets of traits put in are the same, so I shall proceed with the mega models and verify with Lauren
trait_vectors <- as.data.frame(vectors.nmds.comp.trait$vectors$arrows)
NMDS model: only Height, Leaf area, and RMF () found significant.
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = nut_trt, color = nut_trt)) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
# stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Nutrient Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
) +
geom_segment(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2),
inherit.aes = FALSE, # 🚀 Prevents it from looking for 'nut_trt'
arrow = arrow(length = unit(0.2, "cm")),
color = "red",
size = 1
) +
geom_text(
data = as.data.frame(vectors.nmds.comp.trait$vectors$arrows),
aes(x = NMDS1, y = NMDS2, label = rownames(vectors.nmds.comp.trait$vectors$arrows)),
inherit.aes = FALSE, # Avoids unwanted inheritance
hjust = -0.2,
vjust = -0.2,
color = "red",
size = 5
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Year")
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Precipitation Treatment") +
custom_colors +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
geom_point(shape = 21, size = 4, stroke = 0.5) +
stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
labs(title = "Grouped by Grazing History")
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(grazing_hist), fill = as.factor(yr))) +
geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(yr)),
geom = "polygon",
fill = NA,
size = 1) +
labs(title = "Grouped by Grazing History & Year") +
scale_color_manual(
values = c(
"low" = "cyan",
"high" = "brown",
"2019" = "salmon",
"2020" = "limegreen",
"2021" = "cornflowerblue")
)
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
geom_point(shape = 21, size = 4, stroke = 1) +
stat_ellipse(aes(color = as.factor(nut_trt)),
geom = "polygon",
fill = NA, # Makes the ellipse transparent inside
size = 1) +
labs(title = "Grouped by Nutrient Treatment & Year") +
scale_color_manual(
values = c(
"control" = rgb(220, 220, 220, maxColorValue = 255), # Lighter Gray
"+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
"+compost" = rgb(220, 180, 60, maxColorValue = 255), # Bright Golden Yellow
"wet" = rgb(0, 0, 255, maxColorValue = 255),
"dry" = rgb(183, 65, 14, maxColorValue = 255)
)
)